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Many features of granular media can be modelled as a fluid of hard spheres 
with inelastic collisions. Under rapid flow conditions, the macroscopic be- 
havior of grains can be described through hydrodynamic equations. At low- 
density, a fundamental basis for the derivation of the hydrodynamic equations 
and explicit expressions for the transport coefficients appearing in them is pro- 
vided by the Boltzmann kinetic theory conveniently modified to account for 
inelastic binary collisions. The goal of this chapter is to give an overview of the 
recent advances made for binary granular gases by using kinetic theory tools. 
Some of the results presented here cover aspects such as transport properties, 
energy noncquipartition, instabilities, segregation or mixing, non-Newtonian 
behavior, .... In addition, comparison of the analytical results with those ob- 
tained from Monte Carlo and molecular dynamics simulations is also carried 
out, showing the reliability of kinetic theory to describe granular flows even 
for strong dissipation. 

1 Introduction 

Granular systems have attracted the attention of the physics community in 
the past few years, in part because the behavior of these systems under many 
conditions exhibit a great similarity to ordinary fluids pQ. These conditions 
include rapid, dilute flows where the dominant transfer of momentum and 
energy is through binary collisions of the grains. The main difference from 
ordinary fluids is the absence of energy conservation, leading to both obvious 
and subtle modifications of the usual macroscopic balance equations as well as 
the constitutive equations for the irreversible fluxes. However, in spite of the 
utility of hydrodynamics to describe rapid granular flows, there are still some 
open questions about its domain of validity and the associated constitutive 
equations appearing in the hydrodynamic equations [2]. 

To isolate the effects of such collisional dissipation from other important 
properties of granular media, an idealized microscopic model system is usually 
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considered: a system composed by smooth hard spheres with inelastic colli- 
sions. As in the elastic case, the collisions are specified in terms of the change 
in relative velocity at contact but with a decrease in the magnitude of the nor- 
mal component measured by a positive coefficient of restitution a < 1. This 
parameter distinguishes the granular fluid (a < 1) from the ordinary fluid 
(a = 1). Given that the hard sphere system with elastic collisions has been 
widely studied for both equilibrium and non-equilibrium statistical mechanics 
[3], it is tempting to apply the same methods for the case of inelastic collisions. 
However, some care is warranted in translating properties of ordinary fluids to 
granular fluids. In this presentation, a kinetic theory description based on the 
Boltzmann kinetic equation (which applies at sufficiently low density) will be 
considered as the appropriate tool to study granular flows from a microscopic 
point of view. 

Although many efforts have been devoted in the past few years to the 
understanding of granular fluids, the derivation of the form of the constitu- 
tive equations with explicit expressions for the transport coefficients is still a 
subject of interest and controversy. The conditions to obtain a hydrodynamic 
description are expected to be similar to those for normal fluids. For a given 
initial state there are two stages of evolution. First, during the kinetic stage 
there is rapid velocity relaxation to a "universal" velocity distribution that 
depends on the average local density, temperature, and flow velocity. Subse- 
quently, the hydrodynamic stage is described through a slower evolution of 
these local hydrodynamic fields as they approach uniformity. The solution to 
the Boltzmann equation in this second stage is said to be normal, where all 
space and time dependence of the distribution function occurs through the 
macroscopic hydrodynamic fields. The Chapman-Enskog method [4] provides 
a constructive means to obtain an approximation to such a solution for states 
whose spatial gradients are not too large. In this case, the explicit form of this 
normal solution is given as a perturbation expansion in the spatial gradients of 
the fields. This solution is then used to evaluate the fluxes in the macroscopic 
balance equations in terms of these gradients. To lowest order the balance 
equations become the granular Euler equations while to second order they are 
the granular Navier-Stokes equations. In carrying out this analysis, explicit 
forms for the transport coefficients are obtained as functions of the coefficient 
of restitution and other parameters of the collision operator. In this general 
context, the study of hydrodynamics for granular gases is the same as that 
for ordinary fluids. 

The derivation of hydrodynamics from the inelastic Boltzmann equation 
has been widely covered in the case of a monodisperse gas where the particles 
are of the same mass and size. As for elastic collisions, the transport coeffi- 
cients are given in terms of the solutions of linear integral equations [5J [5] , 
which are approximately solved by using Sonine polynomial expansions. The 
estimates for the transport coefficients provided by the Sonine solution com- 
pare in general quite well with both direct Monte Carlo simulation (DSMC) 
of the Boltzmann equation and molecular dynamics (MD) simulation of the 
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gas, even for relatively strong degrees of dissipation [7J [51 M HH]- This good 
agreement supports the formal theoretical analysis and the claim that hydro- 
dynamics is not limited to nearly elastic particles [TTJ [TJl [T3] . 

Nevertheless, a real granular system is generally characterized by some 
degrees of polydispersity in density and size, which leads to phenomena very 
often observed in nature and experiments, such as separation or segregation. 
Needless to say, the study of granular mixtures is much more complicated 
than for a monodisperse gas since not only the number of transport coeffi- 
cients in a multicomponcnt system is higher than that of a single gas, but 
also they depend on parameters such as masses, sizes, composition as well 
as several independent coefficients of restitution a^. Due to these difficulties, 
studies for multicomponcnt gases are more scarce in the literature. Many of 
the previous attempts [14] to derive hydrodynamics from kinetic theory were 
carried out in the quasi-elastic limit where the equipartition of energy can be 
considered as an acceptable assumption. In addition, according to this level of 
approximation, the inelasticity is only accounted for by the presence of a sink 
term in the energy balance equation, so that the expressions for the transport 
coefficients are the same as those obtained for ordinary fluids. However, the 
theoretical prediction of the failure of energy equipartition in multicomponent 
granular gases [15] has been confirmed by computer simulations [16] , and even 
observed in real experiments [17j . 

Although the possibility of nonequipartition was already pointed out many 
years ago [18] , it has not been until recently that a systematic study of the 
effect of nonequipartition on the Navier-Stokes hydrodynamic equations has 
been carefuly carried out 19, 20 . These new equations and associated trans- 
port coefficients provide a somewhat more stringent test of the analysis since 
the parameter space is much larger. As in the monodisperse case, explicit ex- 
pressions for the transport coefficients requires also to consider Sonine poly- 
nomial expansions. The numerical accuracy of this Sonine expansion has been 
confirmed by comparison with Monte Carlo simulations of the Boltzmann 
equation in the cases of the shear viscosity [H] and the tracer diffusion [2"2"] 
coefficients. Exceptions to this agreement are extreme mass or size ratios and 
strong dissipation, although these discrepancies between theory and simula- 
tion diminish as one considers more terms in the Sonine polynomial approxi- 
mation [22] . 

The explicit knowledge of the Navier-Stokes transport coefficients allows 
quantitative application of the nonlinear hydrodynamic equations to a num- 
ber of interesting problems for granular mixtures, such as to quantify the 
violation of the Einstein relation [23, 24 or the Onsager reciprocal relations 
[20] , the stability analysis of the homogeneous cooling state [20] , and segrega- 
tion induced by a thermal gradient [5S]. In all the cases, the analysis clearly 
shows the important role played by the inelasticity in the different physical 
situations. 

The analogy between rapid granular flow and ordinary fluids can be also 
extended to many other transport situations. A particularly simple case, al- 
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lowing detailed analysis even in far from equilibrium conditions is the simple 
or uniform shear flow (USF) problem. Macroscopically, it is characterized by 
uniform density and temperature and a constant mean velocity profile. This 
is a well-known non-equilibrium problem widely studied, for both granular 
monodisperse jT3j [26l [27] and ordinary gases [28]. Nevertheless, the nature of 
this state is quite different in each system. While for elastic fluids the temper- 
ature increases monotonically in time due to viscous heating, a steady state 
is possible for granular media when the effect of the viscous heating is exactly 
compensated by the dissipation in collisions. Thus, in the steady state, there 
is an intrinsic connection between the shear field and dissipation so that the 
collisional cooling sets the strength of the velocity gradient. As a consequence, 
the USF state is inherently non-Newtonian and the rheological properties of 
the system cannot be obtained from the Navier-Stokes description, at least 
for finite dissipation [2T)j . 

The aim of this chapter is to offer a short review of recent results ob- 
tained for binary granular mixtures from the Boltzmann kinetic theory. It is 
structured as follows. The Boltzmann kinetic equation for a granular binary 
mixture and its associated macroscopic balance equations are introduced in 
Sec. [2] Section [3] deals with the solutions to the Boltzmann equation for homo- 
geneous states in the free cooling case as well as when the mixture is heated 
by an external thermostat. The Chapman-Enskog method around the local 
version of the homogenous distributions obtained in Sec. [2] is applied in Sec. [3] 
to get the form of the Navier-Stokes hydrodynamic equations. Theoretical re- 
sults for the diffusion and shear viscosity transport coefficients are compared 
with simulation data in Sec. while the Einstein and the Onsager relations for 
granular mixtures are analyzed in Sees. [6] and [3 respectively. The dispersion 
relations for the hydrodynamic equations linearized about the homogeneous 
cooling state are obtained in Sec. EJ showing that the homogeneous reference 
state is unstable to long wavelength perturbations. The conditions for sta- 
bility are identified as functions of the wave vector, the dissipation, and the 
parameters of the mixture. Segregation due to thermal diffusion is studied in 
Sec. [9] by using the Navier-Stokes description. A new criterion for segregation 
is found that is consistent with recent experimental results. Section fTOl deals 
with the USF problem for a granular mixture. Finally, the paper is closed in 
Sec. [11] with a discussion of the results presented here. 

Before ending this section, I want to remark that the present account is a 
personal perspective based on the author's work and that of his collaborators 
so that no attempt is made to include the extensive related work of many 
others in this field. The references given are selective and apologies are of- 
fered at the outset to the many other important contributions not recognized 
explicitly. 
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2 Boltzmann Kinetic Equation for Binary Mixtures of 
Inelastic Hard Spheres 

Consider a binary granular mixture composed by smooth inelastic disks 
(d = 2) or spheres [d = 3) of masses m\ and m 2 , and diameters o~\ and (7 2 - The 
inelasticity of collisions among all pairs is characterized by three independent 
constant coefficients of restitution an, a 22 , and ai 2 = a 2 i, where a^- < 1 is 
the coefficient of restitution for collisions between particles of species i and j. 
Since the spheres are assumed to be perfectly smooth, only the translational 
degrees of freedom of grains are affected by dissipation. In the low density 
regime, a simultaneous interaction of more than two particles is highly unlike 
and so can be neglected. Consequently, in a dilute gas the interactions among 
the particles reduce to a succession of binary collisions. At this level of de- 
scription, all the relevant information on the state of the system is contained 
in the one-body velocity distribution functions /j(r, v;£) (i = 1,2) defined so 
that fi(r, v; i)drdv is the most probable (or average) number of particles of 
species i which at time t lie in the volume element dr centered at the point 
r and moving with velocities in the range dv about v. For an inelastic gas, 
the distributions /j(r,v;t) (i — 1,2) for the two species satisfy the coupled 
nonlinear Boltzmann equations [301 131] 

U + v • V + ^ • ^ + ^ Mr, v,t)=J2 Jij [v|/i(t), fM , (1) 

where the Boltzmann collision operator [v|/j, fj] is 

Jij [vi|/i, fj] = J dv 2 J dS- 9(a ■ gi 2 )(£ • vu) 

x[ar.Vi(r ) vi,i)/ i (r,v / a ,t)-/i(r,v lj i)/ i (r,v a ,t)] . (2) 

In Eq. ([2]), d is the dimensionality of the system, cr^ = (<jj + aj) /2, a is an 
unit vector along the line of centers, is the Heaviside step function, and 
v i2 = v i — v 2 is the relative velocity. The primes on the velocities denote the 
initial values {v'^v^} that lead to {v 1 ,v 2 } following a binary (restituting) 
collision: 

vi = Vi - fiji (l + a' 1 ) (er • vi 2 )er, 

v 2 =v 2 +^„- (l + a^ 1 ) (ct-Vi2)ct, (3) 

where fiij = mi/ (m, -1- raj). In addition, Fj denotes an external conservative 
force acting on species i (such as a gravity field) and Ti is an operator repre- 
senting a possible effect of an external nonconservative forcing which injects 
energy into the system to compensate for the energy dissipated by collisional 
cooling. This type of force acts as a thermostat that tries to mimics a thermal 
bath. Some explicit forms for the operator T% will be chosen later. 
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The relevant hydro-dynamic fields for the mixture are the number densities 
Ui, the flow velocity u, and the temperature T. They are defined in terms of 
moments of the velocity distribution functions fi as 



dv/i(v) , (4) 

pu = X! m ' / dvv /»( v ) > ( 5 ) 

1=1 ^ 

nT = V = Y j n l T l = Y j ^ |dvF 2 /i(v) , (6) 

z=l z=l 

where V = v u is the peculiar velocity, n = n\ + 71.2 is the total number 
density, p — mini + 7712^2 is the total mass density, and p is the pressure. 
Furthermore, the third equality of Eq. ^ defines the kinetic temperatures Tj 
for each species, which measure their mean kinetic energies. 

The collision operators conserve the particle number of each species and 
the total momentum, but the total energy is not conserved: 

JdvJ iJ [v\f i ,f j ] = 0, (7) 

rrn J dwJy[v|/i, fj] + m 3 J dwJji[v\fj, fi] = , (8) 
22 

Y,Y, m - dv ^% ma, fi) = ~ dnT c • ( 9 ) 

„- — 1 A — 1 J 



i=l j=l 



In Eq. ([9|), (, is identified as the total "cooling rate" due to collisions among 
all species. It measures the rate of energy loss due to dissipation. At a kinetic 
level, it is also convenient to introduce the "cooling rates" Q for the partial 
temperatures Ti. They are defined as 



2 2 



C» = = J dWVtfM/,,^], (10) 

where the second equality defines the quantities Qj ■ The total cooling rate £ 
can be written in terms of the partial cooling rates Q a s 

2 

C = T- 1 ^x i T i G, (11) 

i=l 

where x% — rii/n is the mole fraction of species i. 

From Eqs. ([T])-©, the macroscopic balance equations for the number den- 
sities rij, the total momentum density pxi and the energy density (d/2)nT can 
be obtained. They are given, respectively, by |19j 
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D t m + n,V ■ u + = , (12) 

m l 

1 

^tU + V-P = ^n ! F ! , (13) 

i=l 

DtT — + #" (Vq+PiVu-^T^^U -(C-OT. (14) 
In the above equations, D t — dt + u • V is the material derivative, 

ji =mi J dvV/,(v) (15) 
is the mass flux for species i relative to the local flow, 

P = J2 m * / dWV/ 4 (v) (16) 
is the total pressure tensor, and 

q = ET/ dv ^ 2v/i(v) (17) 

i=l 

is the total heat flux. On the right-hand side of the temperature equation 
([HI) , the source term £ (measuring the rate of heating due to the external 
thermostat) is given by 

1 2 r 

t = -d^T S m V dv ^ 2 ^/i( v )' ( 18 ) 

2=1 

In the balance equations (fI"2"l) - (fT4")) it is assumed that the external driving 
thermostat does not change the number of particles of each species or the 
total momentum, i.e., 

J dv^/itv) = 0, (19) 
Y. m * j dvv^/i(v) =0. (20) 



The utility of the balance equations (JT2J) — (JT4J) is limited without further 
specification of the fluxes and the cooling rate, which in general have a com- 
plex dependence on space and time. However, for sufficiently large space and 
time scales, one expects that the system reaches a hydrodynamic regime in 
which all the space and time dependence is given entirely through a functional 
dependence on the six hydrodynamic fields n.j, u, and T. The corresponding 
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functional dependence of j^, P, and q on these fields are called constitutive 
equations and define the transport coefficients of the mixture. The primary 
feature of a hydrodynamic description is the reduction of the description from 
many microscopic degrees of freedom to a set of equations involving only six 
local fields. At a kinetic level, the constitutive equations are obtained when 
one admits the existence of a normal solution to the Boltzmann equation 
where the velocity distribution functions depend on r and t only through 
their functional dependence on the fields, namely, 

/i(r,vi,*) = Mvi\n t (r,t),T(r,t),u(r,t)}. (21) 

This normal solution is generated by the Chapman-Enskog method [4] con- 
veniently adapted to dissipative dynamics. Since the method is based on an 
expansion around the local version of the homogeneous state, let us charac- 
terize it before considering inhomogeneous solutions. 

3 Homogeneous States 

In this Section we are interested in spatially homogeneous isotropic states. In 
this case, we assume that the magnitude of the conservative external forces 
is at least of first order in the spatial gradients (i.e., F< = 0), so that Eq. ([I]) 
becomes 

(d t +Fi)fi=J2 MfiJjl (22) 

3 

For elastic collisions (a^ = 1) and in the absence of external forcing {Ti = 0), 
it is well known that the long-time solution to ([22]) is a Maxwellian distribution 
for each species at the same temperature T. However, if the particles collide 
inelastically {otij < 1) and T = 0, a steady state is not possible in uniform 
situations since the temperature decreases monotonically in time. This state is 
usually referred to as the homogeneous cooling state (HCS). In this situation, 
since rii is uniform and u = 0, the normal (hydrodynamic) solution to /, 
requires that all its time dependence occurs only through the temperature 
T(t). Consequently, fi(v,t) must be of the form [T5] 

/i(v,t) =mvo d (m(v/v (t)) , (23) 

where vq (t) = y / 2T(t)(mi + rnv)/ (m\m2) is a thermal speed defined in terms 
of the temperature T(t) of the mixture. The balance equations (fT2")) - (ri"4")) to 
this order become dtXi = dtu = 0, and T _1 9 t T = — £, where the cooling rate 
C is determined by Eq. ([9]). In addition, from Eqs. (fT0|) and (l22l) (when Ti = 0) 
one can derive the time evolution for the temperature ratio 7 = Ti(t)/T2(t): 



9t In 7 = £2 - Ci • 



(24) 
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The fact that the distributions depend on time only through T(t) necessarily 
implies that the temperature ratio 7 must be independent of time and so, Eq. 
(pM|) gives the condition 

(i=(2= C (25) 

In the elastic case, where fi is a Maxwellian distribution, the above condition 
yields Ti(i) = T2(t) = T(t) and the energy equipartition applies. However, in 
the inelastic case, the equality of the cooling rates leads to different values for 
the partial temperatures, even if one considers the Maxwellian approximation 
to fi. Nevertheless, the constancy of 7 assures that the time dependence of 
the distributions is entirely through T(t), and in fact the partial temperatures 
can be expressed in terms of the global temperature as 

r ^= l + J 7 -D rW ' ^ l + ai ( 7 -l) T(t) - (26) 



Just as for the single gas case [30 , 32 , the exact form of <Pi has not yet been 
found, although a good approximation for thermal velocities can be obtained 
from an expansion in Sonine polynomials |15j . In the leading order, <Pi is given 

by 



Here, v* = v/vq, 



(27) 



2 

m 



V 3=1 

and 7^ = Ti/T. The coefficients Cj (which measure the deviation of <Pi from 
the reference Maxwellian) are determined consistently from the Boltzmann 
equation. The explicit form of the approximation (|27p provides detailed pre- 
dictions for the temperature ratio T1/T2 (through calculation of the cooling 
rates) and for the cumulants c$ as functions of the mass ratio, size ratio, com- 
position and coefficients of restitution |15j . The numerical accuracy of this 
truncated Sonine expansion has been confirmed by comparison with Monte 
Carlo [33] and MD [33] simulations. 

As said in the Introduction, the existence of different temperatures for 
each species has been observed in real experiments of driven steady states. 
These states are achieved from external forces that do work at the same rate 
as collisional cooling. In experiments this is accomplished by vibrating the 
system so that it is locally driven at walls. Far from these walls a steady state 
is reached whose properties are expected to be insensitive to the details of 
the driving forces. Due to the technical difficulties involved in incorporating 
oscillating boundary conditions, it is usual to introduce external forces (or 
thermostats) acting locally on each particle. These forces are represented by 
the operator T\ in Eq. (|22p and depend on the state of the system. Two 
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types of thermostats have been usually considered in the literature. One is a 
deterministic thermostat widely used in nonequilibrium MD simulations for 
ordinary fluids [35] [M]- The force is similar to a Stokes law drag force, linear 
in the velocity, but with the opposite sign so that it heats rather than cools 
the system. In this case, Ti is given by [2TJ [37] 



Fifi(v) = ■ [v/i(v)], (29) 

where the friction constant in the force has been adjusted to get a constant 
temperature in the long-time limit. It must be remarked that the correspond- 
ing Boltzmann equation (|22[) for this Gaussian thermostat force is formally 
identical with the Boltzmann equation in the HCS (i.e., with T% = 0) when 
both equations are written in terms of the reduced distributions <Pi(v*) [37j . 
In particular, the dependence of 7 on the parameters of the system is the same 
with and without the Gaussian thermostat. 

A second method of driving the system is by means of a stochastic 
Langevin force representing Gaussian white noise |38j . The corresponding op- 
erator Ti has a Fokker-Planck form [32] 

IT, / d x 2 



^(v)=-~C^j *(v), (30) 

where for simplicity the covariance of the stochastic acceleration has been 
taken to the same for each species [39[ [40] . This requirement gives the steady 
state condition 

Ti To 

— Ci = — C2. (31) 
mi m2 

The cooling rates Q are no longer equal, in contrast to the HCS, and the de- 
pendence of 7 on the control parameters is different as well [34] . The procedure 
for determining the temperature ratio and the cumulants Ci is the same as in 
the HCS state since the steady state distribution <Pi can also be represented as 
an expansion of the form (|27p and the coefficients are now determined from 
the solution to the Boltzmann equation ([2"2"]) . The condition (13T]) gives the 
corresponding equation for the temperature ratio. 

Figure [T] illustrates the differences between the HCS and the stochastic 
steady state at the level of the temperature ratio T1/T2. We have considered 
mixtures constituted by spheres (d = 3) of the same material [and so, otij = a, 
and rai/7712 — (C1/CT2) 3 ] and equal volumes of large and small particles [i.e., 
%2 = (ci/o^) 3 ^]- Here, for the sake of simplicity, the cooling rates have been 
analytically estimated by using Maxwellians (namely, by taking Cj = 0) for 
the distributions /t(v): 



in (d-l)/2 2 /a , fl N i /2 



x(l + ay) 



(32) 
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Fig. 1. Temperature ratio T\jT% versus the size ratio 0-1/02 for a.ij = a — 0.78 in 
the case of mixtures constituted by particles of the same material and equal total 
volumes of large and small particles. The lines are the kinetic theory results in (a) 
the stochastic driving case and (b) the free cooling case, while the points refer to 
MD simulations [41] . 

Simulation data recently obtained from MD simulations in agitated mixtures 
have also been included [41] . The experimental value of the (common) coef- 
ficient of restitution is a = 0.78. While a good agreement between kinetic 
theory and MD simulations is found when the gas is assumed to be driven 
by the stochastic thermostat, significant discrepancies appear in the undriven 
(HCS) case, especially as the size ratio o"i/<72 increases. These results con- 
trast with the comparison made by Brey et al. [42] for agitated mixtures in 
the tracer limit (x\ — > 0), where the predictions of the temperature ratio from 
kinetic theory based on the condition £i = £ 2 compare quite well with MD 
simulations. However, for the cases studied in Ref. [H] the conditions (p?5]) and 
(pTT]) yield quite similar results for the dependence of Tx/T 2 on the parameters 
of the system. The good agreement found in Fig. [1] between MD simulations 
for agitated mixtures and kinetic theory suggests that the stochastic driving 
condition can be considered as a plausible first approximation for qualitative 
comparisons with experimental results [17] . 

4 Navier— Stokes Hydrodynamic Equations 

We consider now a spatially inhomogeneous state created either by initial 
preparation or by boundary conditions. We are interested in a hydrodynamic 
description where the state of the system is completely specified through their 
hydrodynamic fields. This implies that the latter dominate over other excita- 



12 Vicente Garzo 



tions for times large compared to the mean free time and for wavelengths large 
compared to the mean free time. The hydrodynamic regime is characterized 
by the existence of a normal solution to the Boltzmann equation which can be 
explicitly obtained by means of the Chapman-Enskog method [?]. For small 
spatial variations, the functional dependence (f2Tj) of the normal solution can 
be made local in space and time through an expansion in gradients of the 
fields: 

iWf 0) + e/f ) + ---, (33) 

where each factor of e means an implicit gradient of a hydrodynamic field. 
The reference distribution function (r,v, i) is the local version of the ho- 
mogeneous distribution (|23[) . namely, it is obtained from the homogeneous 
distribution by replacing the temperature, densities and flow velocity by their 
noncquilibrium values: 

/i 0) (r,v,t) = n i (r,f) Uo - d (T(r,t))^ (V/v (T(r,t))) , (34) 

where V = v — u(r, t). The time derivatives of the fields are also expanded as 
d t = + edf + • • • . The coefficients of the time derivative expansion are 
identified from the balance equations fT/T|) - (fr4")) with a representation of the 
fluxes and the cooling rate in the macroscopic balance equations as a similar 
series through their definitions as functionals of 

This is the usual Chapman-Enkog method for solving kinetic equations 
[HE3H3]. Nevertheless, the complexity introduced by the energy dissipation 
in collisions has led to the introduction by some authors [H) of some additional 
approximations, restricting the validity of most of the results to the small 
inelasticity limit. Only very recently, explicit expressions for the fluxes to first 
order in the gradients as explicit functions of the coefficients of restitution have 
been obtained [HI EH [44] . To Navier-Stokes order, the constitutive equations 
for mass, momentum, and heat fluxes are given, respectively, by 



2 




2 

q (i) = _ T 2 D " Vxi _ LVp _ XVT + KiFi. (37) 

The transport coefficients in these equations are the diffusion coefficient D, 
the pressure diffusion coefficient D p , the thermal diffusion coefficient D' , the 
mobility coefficient Xiji the- shear viscosity rj, the Dufour coefficient D", the 
pressure energy coefficient L, the thermal conductivity A, and the coefficient 
Hi. These coefficients are defined as 
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J dvV- A u 



^/.IvV.B, 

L y = /dvv-d, 



7- 



\ ^ / dvV-£ ij; 



(d-l)(d+2)£- 



1 2 
dT 2 ^ 



i=l 



dvVV : T> u 



dvV 2 V ■ A, 



dv V 2 V ■ 



X = 



i=\ J 

-t 



2 
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(38) 

(39) 
(40) 
(41) 

(42) 
(43) 
(44) 
(45) 
(46) 



Here, .Aj(V), Bj(V), Cj(V), 2?j(V), and £ij(V) are functions of the peculiar 
velocity and the hydrodynamic fields. They obey the following set of coupled 
linear integral equations: 



4 0) {Td T + pd p ) + F\ u> + & Ai + M X A 2 = Ax 



r(0) 



dc 



(0)' 



dxi 



x(pSi+TCi), (47a) 



4 0) (Ta T +pa p ) + ^ UJ + /:2 



(0) 



A 2 + M 2 A X = A 2 + I 



x( P 5 2 + TC 2 ), (47b) 



c[ 0) (Ta T +p9 P )+^r+A 



(o) 



2c 



(0) 



h+MiB 2 = Bi 



Tc 



(o) 



-Ci, (48a) 



c [ 0) (Ta T +pa p ) + ^ lv +/: 2 -2c^ b 2 +m 2 Bi = b 2 +^^c 2 , (48b) 



-(0) 



JO) 



Tc 



(0) 
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cf ] (Td T + pd p ) + F[ Q) + d 



-4 0) 

2 c 



c 1 +m 1 c 2 = c 1 



pel 



(0) 



2T 



■Si, (49a) 



(0) 



-c (0) 



(o) 

C 2 + M 2 d =03-^62, (49b) 



-4 Uj (Ta T + pa p ) + ^ 0) + /: 1 
-4 0) (Ta T + P a p ) + ^ 0) + /: 2 



-4 0) (TS T + p9 p )+^ 0) +A 
-4 0) (T8t + pdp) + jf> + & 



T> 1 +MiD 2 = D 1 , (50a) 

X» 2 + M 2 T> 1 = D 2 , (50b) 

fii+Mif 2 i = E n , (51a) 

£i2+Xi5 22 =Ei 2 . (51b) 



Here, we have introduced the linearized Boltzmann collision operators 
CiX = - (jnlA°\x] + Ju[X, / x (0) ] + J a [X, /f] 

MiX = -J 12 [ri°\x], 

with a similar definition for C 2 and Ai 2 . In addition, = in the undriven 

case while = in the driven case, where is given by Eq. |9]) to zeroth 
order, i.e., 



(52) 



(53) 



(0) 



^EE m ,/dvy%-[v|f,/f] 



»=i j=i 

In Eqs. I|47ap - (j51b[) we have also introduced the operators 



TX''X = -Cr J ^ • (VI) , (Gaussian thermostat), 

Tf^X = — — — Ci°' ) ( J (stochastic thermostat), 

and the quantities 



Ai(V) 



_a_ f ( ) 



v. 



Bi(V) 



Ci(V) 



nT ( d 



f(0) 



T 



7i 2 9V 



(v/f) 



v. 



(54) 

(55a) 
(55b) 

(56) 
(57) 
(58) 
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D i (v)=vA / (o) -ifv-A/CoA,, (59) 



dV l d\ dV 



W = - s?f rK-? ' (60) 



Here, I is the unit tensor in d dimensions and pi — mini is the mass density 
of species i. Upon writing Eqs. (|47a|) — (|51b[) use has been made of the fact 
that there is no contribution to the cooling rate at this order, i.e., £W = 0. 
As a consequence, J-^ — 0. The property = is special of the low 
density Boltzmann kinetic theory (since /| does not contain any contribution 
proportional to V • u) , but such terms occur at higher densities [151 135] . Note 
that in the particular case of the gravitational force Pj = m^g (where g is 
the gravity acceleration), the combination miEn +m2Ei2 = 0. This leads to 
£ij = 0, and so there are no contributions to the mass and heat fluxes coming 
from the external conservative forces. 

In summary, the solutions to the Boltzmann equations to first order in the 
spatial gradients are given by [19j 

fi = /■ ° } +A i -Vx 1 + B i -Vp + C i -VT + T>i :Vu + ^ S l3 -F r (61) 



The solution to zeroth-order is obtained from Eq. (|22p while the functions 
{A-i, Bi, Ci, T>i, £ij} are determined from the integral equations (|47a|) — (|51b|l . 
Once these equations are solved, the Navier-Stokes transport coefficients are 
obtained from Eqs. (|58)) - (|4"6"]) and the mass, momentum, and heat fluxes are 
explicitly known. These fluxes, together with the macroscopic balance equa- 
tions (fT2"]) - (fT4")) . provide the closed set of Navier-Stokes order hydrodynamic 
equations for a granular binary mixture. All these results are still formally 
exact and valid for arbitrary values of the coefficients of restitution. 

However, explicit expressions for the Navier-Stokes transport coefficients 
require to solve the integral equations (|47a|) — (|51b|l . Accurate approximations 
for {>ti, Bi, Ci, 2?i, £ i 3 } may be obtained using low order truncation of ex- 
pansions in a series of Sonine polynomials. The polynomials are defined with 
respect to a Gaussian weight factor whose parameters are chosen such that 
the leading term in the expansion yields the exact moments of the entire dis- 
tribution with respect to 1, v, and v 2 . The procedure is similar to the one 
followed for elastic collisions [1] and yields explicit expressions for the trans- 
port coefficients in terms of the parameters of the mixture [TH1 HOI El] • 



5 Comparison with Monte Carlo Simulations 

As said before, the expressions derived for the Navier Stokes transport coef- 
ficients are obtained by considering two different approximations. First, since 
the deviation of ff 3 from its Maxwellian form is quite small in the region of 
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thermal velocities, one uses the distribution l|27p as a trial function for ' . 
Second, one only considers the leading terms of an expansion of the distri- 
bution f^ 1 ' in Sonine polynomials. Both approximations allow one to offer 
a simplified kinetic theory for a granular binary mixture. To assess the de- 
gree of accuracy of these predictions, one resorts to numerical solutions of the 
Boltzmann equation, such as those obtained from the Direct simulation Monte 
Carlo (DSMC) method [17]. Although the method was originally devised for 
normal fluids, its extension to granular gases is straightforward [37] . In this 
Section we provide some comparisons between theory and numerical solutions 
of the Boltzmann equation by means of the DSMC method in the cases of the 
diffusion coefficient D (in the tracer limit) and the shear viscosity coefficient 
rj of a heated gas. Let us study each coefficient separately. 

5.1 Tracer Diffusion Coefficient 

We consider a free granular mixture {J~i = 0) in which one of the components 
of the mixture (say, for instance, species 1) is present in tracer concentration 
(xi — > 0). In this situation the diffusion coefficient of impurities in a granular 
gas undergoing homogeneous cooling state can be measured in simulations 
from the mean square displacement of the tracer particle after a time interval 
t 08]: 

— <|r(t) - r(0)| 2 ) = , (62) 

at n>2 

where |r(t) — r (0) | is the distance travelled by the impurity from t = until 
time t. The relation (|62[) written in appropriate dimensionless variables to 
eliminate the time dependence of D(t) can be used to measure by computer 
simulations the diffusion coefficient [HI [22 • 

If the hydrodynamic description (or normal solution in the context of the 
Chapman-Enskog method) applies, then the diffusion coefficient D(t) depends 
on time only through its dependence on the temperature T(t). Dimensional 
analysis shows that D(t) oc y/T(t). In this case, after a transient regime, 
the reduced diffusion coefficient D* = {mim^J 'p)D(t)vo(t)/T(t) achieves a 
time-independent value. Here, vo(t) = naf^^oit) oc y/T(t) is an effective 
collision frequency for hard spheres. The fact that D* reaches a constant 
value for times large compared with the mean free path is closely related 
with the validity of a hydrodynamic description for the system. In addition, 
as has been recently shown [49], the dependence of D* on the mass ratio 
7Tii/m2 and the coefficient of restitution olxi is only through the effective 
mass m* = m\ + (mi + m 2 ){l — Qi2)/(1 + "12)- 

The dependence of D* on the common coefficient of restitution = a 
is shown in Fig. [5] in the case of hard disks (d — 2) for three different sys- 
tems. The symbols refer to DSMC simulations while the lines correspond to 
the kinetic theory results obtained in the first Sonine approximation |19[ I44j . 
MD results reported in Ref. [8] when impurities and particles of the gas are 
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Fig. 2. Plot of the reduced diffusion coefficient D* as a function of the (common) 
coefficient of restitution a for binary mixtures with uj = fi in the case of a two- 
dimensional system (d = 2). Here, ui = <ti/<T2 and /i = mi/m2. The symbols are 
computer simulation results obtained from the mean square displacement and the 
lines are the theoretical results obtained in the first Sonine approximation. The 
DSMC results correspond to ^ = 1/4 (•), /x = 4 (O) and fi = 1 (O). Molecular 
dynamics results reported in Ref. [8] for [i = 1 (A) have also been included. 

mechanically equivalent have also been included. We observe that in the latter 
case MD and DSMC results are consistent among themselves in the range of 
values of a explored. This good agreement gives support to the applicability of 
the inelastic Boltzmann equation beyond the quasielastic limit. It is apparent 
that the agreement between the first Sonine approximation and simulation 
results is excellent when impurities and particles of the gas are mechanically 
equivalent and when impurities are much heavier and/or much larger than 
the particles of the gas (Brownian limit). However, some discrepancies be- 
tween simulation an theory are found with decreasing values of the mass ratio 
fi = mi/m,2 and the size ratio u> = ci/o^- These discrepancies are not easily 
observed in Fig. [5] because of the small magnitude of D* for fi = 1/4. For 
these systems, the second Sonine approximation [22j improves the qualita- 
tive predictions over the first Sonine approximation for the cases in which 
the gas particles are heavier and/or larger than impurities. This means that 
the Sonine polynomial expansion exhibits a slow convergence for sufficiently 
small values of the mass ratio fi and/or the size ratio u. This tendency is also 
present in the case of elastic systems |50j . 
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5.2 Shear Viscosity Coefficient of a Heated Gas 

The shear viscosity rj is perhaps the most widely studied transport coefficient 
in granular fluids. This coefficient can be measured in computer simulations 
in the special hydrodynamic state of uniform shear flow (USF). At a macro- 
scopic level, this state is characterized by constant partial densities rij, uni- 
form temperature T, and a linear flow velocity profile u\ t k = U2.k — akife, 
aki = aSkxSiy, a being the constant shear rate. In this state, the temperature 
changes in time due to the competition between two mechanisms: on the one 
hand, viscous heating and, on the other hand, energy dissipation in collisions. 
In addition, the mass and heat fluxes vanish by symmetry reasons and the 
(uniform) pressure tensor is the only nonzero flux of the problem. The relevant 
balance equation is that for temperature, Eq. (fl"4"|) . which reduces to 

d t T + ^-aP xy = -(£ - OT, (63) 

where 

P xy =Y d m i f dVV x V y fi(V) (64) 
i=i ^ 

is the xy-element of the pressure tensor. 

For a granular fluid under USF and in the absence of a thermostatting 
force (£ = 0), the energy balance equation (|63|) leads to a steady state when 
the viscous heating effect is exactly balanced by the collisional cooling. This 
situation will be analyzed in Sec. 1101 However, if for instance the mixture 
is heated by the Gaussian thermostat (|29|) (with v — > V), then the viscous 
heating still prevails so that the temperature increases in time. In this case, 
the collision frequency vq (t) oc y/T{t) also grows with t and hence the reduced 
shear rate a*(t) = a/i>o(t) (which is the relevant nonequilibrium parameter 
of the problem) monotonically decreases in time. Under these conditions, the 
system asymptotically achieves a regime described by linear hydrodynamics 
and the (reduced) shear viscosity rf = \vo{t) / nT (t)\q(t) can be measured as 

p* 

rf = - lim (65) 

where P* y = P xy /nT. This procedure allows one to identify the shear viscosity 
of a granular mixture excited by the Gaussian external force (|29|) and compare 
it with the predictions given by the Chapman-Enskog method. 

In Fig. [3l we plot the ratio rj* (a) / n* (1) versus the mass ratio mi/rri2 in the 
case of hard spheres (d = 3) for (71/02 = 1, X\ = i, and three different values 
of the (common) coefficient of restitution = a. Here, rf(l) refers to the 
elastic value of the shear viscosity coefficient. Again, the symbols represent 
the simulation data obtained by numerically solving the Boltzmann equation 
[21] . while the lines refer to the theoretical results obtained from the Boltz- 
mann equation in the first Sonine approximation. We see that in general the 
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Fig. 3. Plot of the ratio rf (a) / 'rf (1) as a function of the mass ratio m\/m,2 for 
<yi/<J2 = n\/n2 = 1 and three different values of the (common) coefficient of resti- 
tution a: a = 0.9 (circles), a = 0.8 (squares), and a = 0.7 (triangles). The lines are 
the theoretical predictions and the symbols refer to the results obtained from the 
DSMC method. 
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Fig. 4. Plot of the ratio rj* (a) / rf (1) as a function of the size ratio ai/02 for 
ra\jm,2 = 4, ni/n2 = 1 and three different values of the (common) coefficient of 
restitution a: a — 0.9 (circles), a — 0.8 (squares), and a = 0.7 (triangles). The lines 
are the theoretical predictions and the symbols refer to the results obtained from 
the DSMC method. 
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deviation of r/* (a) from its functional form for elastic collisions is quite impor- 
tant. This tendency becomes more significant as the mass disparity increases. 
The agreement between the first Sonine approximation and simulation is seen 
to be in general excellent. This agreement is similar to the one previously 
found in the monocomponent case [7] [TO] [51] . At a quantitative level, the dis- 
crepancies between theory and simulation tend to increase as the coefficient 
of restitution decreases, although these differences are quite small (say, for 
instance, around 2% at a = 0.7 in the disparate mass case m\jm% — 10). 
The influence of the size ratio on the shear viscosity is shown in Fig. 0] for 
mi /ma = 4 and x\ — \ [21]. We observe again a strong dependence of the 
shear viscosity on dissipation. However, for a given value of a, the influence of 
G\ju2 on rf is weaker than the one found before in Fig. [3] for the mass ratio. 
The agreement for both a — 0.9 and a = 0.8 is quite good, except for the 
largest size ratio at a = 0.8. These discrepancies become more significant as 
the dissipation increases, especially for mixtures of particles of very different 
sizes. In summary, according to the comparison carried out in Figs. [3] and [U 
one can conclude that the agreement between theory and simulation extends 
over a wide range values of the coefficient of restitution, indicating the relia- 
bility of the first Sonine approximation for describing granular flows beyond 
the quasielastic limit. 

6 Einstein Relation in Granular Gases 

The results presented in Section [5] give some support to the validity of the 
hydrodynamic description to granular fluids. However, in spite of this sup- 
port some care is warranted in extending properties of normal fluids to those 
with inelastic collisions. Thus, for elastic collisions, in the case of an impurity 
(tracer) particle immersed in a gas the response to an external force on the 
impurity particle leads to a mobility coefficient proportional to the diffusion 
coefficient. This is the usual Einstein relation 48], which is a consequence of 
the fluctuation-dissipation theorem. A natural question is whether the Ein- 
stein relation also applies for granular fluids. 

To analyze it, let us consider the tracer limit [xi — > 0) and assume that 
the current of impurities is only generated by the presence of a weak 
concentration gradient Vxi and/or a weak external field Fi acting only on 
the impurity particles. Under these conditions, Eq. (|3"S"|) becomes 



The Einstein ratio e' between the diffusion coefficient D and the mobility 
coefficient \n is defined as 



jW _ -miDVaji+xuFi- 



(66) 



D 



(67) 



m\X\ 
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where T ~ T 2 111 the tracer limit. For elastic collisions, the Chapman-Enskog 
results yield e' = 1. However, at finite inelasticity the relationship between D 
and xn is no longer simple and, as expected, the Chapman-Enskog expres- 
sions for D and xn i n the case of an unforced granular gas clearly show 
that e' ^ 1. This means that the Einstein relation does not apply in granular 
gases. The deviations of the (standard) Einstein ratio e' from unity has three 
distinct origins: the absence of the Gibbs state (non-Gaussianity of the distri- 
bution function of the HCS) , time evolution of the granular temperature, and 
the occurrence of different kinetic temperatures between the impurity and gas 
particles. The second source of discrepancy can be avoided if the system is 
driven by an external energy input to achieve a stationary state. With respect 
to the third reason of violation, this could also be partially eliminated if the 
temperature of the gas T is replaced by the temperature of the impurity T\ 
in the usual Einstein relation (|67p . This change yields the modified Einstein 
ratio 

D , . 

e = mixi— . (68) 

As a consequence, the only reason for which e ^ 1 is due to the non-Maxwellian 
behavior of the HCS distribution. Given that the deviations of the gas dis- 
tribution from its Maxwellian form are small T5j, the discrepancies of e 
from unity could be difficult to detect in computer simulations. This conclu- 
sion agrees with recent MD simulations [52] of granular mixtures subjected to 
the stochastic driving of the form ([30]) . where no deviations from the (modi- 
fied) Einstein relation e = 1 have been observed for a wide range of values of 
the coefficients of restitution and parameters of the system. 

To illustrate the influence of dissipation on the Einstein ratio more gener- 
ally, in Figs. [5] and [5] the Einstein ratio as given by ([6"5]) is plotted versus the 
coefficient of restitution ai2 for cxx/a^ — 1 and different values of the mass ra- 
tio mi /ni2 and the coefficient of restitution U21 ■ The results obtained by using 
the Gaussian thermostat ([29]) are shown in Fig. [5] while Fig. [6] corresponds 
to the results derived when the system is heated by the stochastic thermostat 
(|5U|) [25]. We observe that in general e / 1, although its value is very close to 
unity, especially in the case of the stochastic thermostat, where the deviations 
from the Einstein relation are smaller than 1%. However, in the case of the 
Gaussian thermostat the deviations from unity are about 8%, which could be 
detected in computer simulations. Figures [5] and [6] also show the fact that the 
transport properties are affected by the thermostat introduced so that the 
latter does not play a neutral role in the problem [51j . 

7 Onsager's Reciprocal Relations in Granular Gases 

In the usual language of the linear irreversible thermodynamics for ordinary 
fluids [53] , the constitutive equations ([35]) and (]37p for the mass flux and heat 
flux in the absence of external forces can be written as 
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Fig. 5. Plot of the modified Einstein ratio e versus the coefficient of restitution cm 
for the Gaussian thermostat in the cases: (a) «22 = cm, mi/m.2 = 5 and o\/a2 = 1; 
(b) «22 = cm, mi/rri2 = 0.5 and cr\/o2 = 1; and (c) 022 = 0.5, mi/m.2 = 10 and 
0-1/0-2 = 1. 
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Fig. 6. Plot of the Einstein ratio e versus the coefficient of restitution cm for the 
stochastic thermostat in the cases: (a) a.22 = cm, mi/rri2 = 5 and ai/02 = 1; (b) 
Q22 = cm, mi/rri2 = 0.5 and ai/a2 = 1; and (c) 022 = 0.5, mi/m.2 = 10 and 
01/0-2 = 1. 
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jf^-E^f^) -^-C P Wp, (69) 



where 



and 



JW = -L M VT-E^(^ 1 ) -ClXl'. (70, 



^) =— Vln(^), (72) 



Hi being the chemical potential per unit mass. In Eqs. (|69[) and ([70]) . the 
coefficients arc the so-called Onsager phenomenological coefficients and the 
coefficients C p and C' p can be expressed in terms of the transport coefficients 
associated with the heat and mass fluxes. For elastic fluids, Onsager showed 
[53] that time reversal invariance of the underlying microscopic equations of 
motion implies important constraints on the above set of transport coefficients, 
namely 

Lij — Lji, Li q — L q i, C'p = C' p = 0. (73) 

The first two symmetries are called reciprocal relations as they relate trans- 
port coefficients for different processes. The last two are statements that the 
pressure gradient does not appear in any of the fluxes, even though it is ad- 
mitted by symmetry. Even for a one component fluid, Onsager's theorem is 
significant as it leads to a new contribution to the heat flux proportional to 
the density gradient [5] . Since there is no time reversal symmetry for granular 
fluids, Eqs. (I73|) cannot be expected to apply. However, since explicit expres- 
sions for all transport coefficients are at hand, the quantitative extent of the 
violation can be explored. 

To make connection with the expressions (j35| and (|37p for the mass and 
heat fluxes, respectively, it is first necessary to transform Eqs. (f69"|) - (|71ll to 
the variables xx, p, T. Since Vxi = — Vx2, Eq. (|72|) implies 



(V/Zi)t - (V^ 2 )t np 



T pip2 



vx± H (m2 — mjJVlnp 

np 



(74) 



The coefficients {L^, L iq , L qi , L qq7 C p ,C' p } then can be easily obtained in 
terms of the Navier-Stokes transport coefficients introduced in Sec. |U The 
result is [HI 

Lxx = -L 12 = — L 2 i = = D, Lx q = pi D , (75) 

P 

T T T 2 Pl p 2 „ d+2 T Pl p 2 , , , > 

Lqi = -L q2 = D — (m 2 - mx)D, (76) 

np 2 p z 
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d + 2 mi - mi , 
2 m\m 2 

C p = ^D p -^(m 2 - mi )D, (78) 

C=L 1 -C p - -±^T 2 m 2 - mi )D" . 79 

' I p m\mi npp 

Onsager's relation Li 2 — L 2 \ holds since the diffusion coefficient D is 
symmetric under the change 1 2 [15]. However, in general L\ q ^ L ql , C p ^ 
0, and C' p ^ PHI- The Chapman- Enskog results [Hj show that there are only 
two limit cases for which L± q — L q \ — C p = C' p = 0: (i) the elastic limit (ay 
1 : with arbitrary values of masses, sizes and composition and (ii) the case 
of mechanically equivalent particles with arbitrary values of the (common) 
coefficient of restitution a = ay. Beyond these limit cases, Onsager's relations 
do not apply. At macroscopic level the origin of this failure is due to the cooling 
of the reference state as well as the occurrence of different kinetic temperatures 
for both species. 



8 Linearized Hydrodynamic Equations and Stability of 
the Homogeneous Cooling State 

As shown in Sec. [4] the Navier-Stokes constitutive equations (|35|) (|37| have 
been expressed in terms of a set of experimentally accessible fields such as 
the composition of species 1, x±, the pressure p, the mean flow field u, and 
the granular temperature T. In terms of these variables and in the absence of 
external forces, the macroscopic balance equations (fT2" |) -lfT^ |) become 

+ ^ V-ji = 0, (80) 

2 

D t p + pV- u+ - (V • q+ P : Vu) = -(p, (81) 
a 

Au + p^V • P = , (82) 

D t T--Y ^ + ^-(V-q+P: Vu) = -(T . (83) 
n mi an 

i 

When the expressions ([35 |) -([37 |l for the fluxes and the cooling rate £ — > Q ^ are 
substituted into the above exact balance equations ([80 |) -([83 |) one gets a closed 
set of hydrodynamic equations for x± , u, T, and p. These are the Navier-Stokes 
hydrodynamic equations for a binary granular mixture: 

D tXl = -J^—^ • ( ^^DV Xl + P -D P V P + ZiyVT) , (84) 
n z mim2 \ p p 1 J 
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(A + c) p + ~^-pV ■ u = 2 v ' ( t2d " Vx i + l ^p + AVT ) 

2 / 2 \ 

(85) 



2 _ T Too — wi „ / mvm^n „„ p „ „ 

A + C T + ^pV • u = -V • DV Xl + -D p S7p 

an n m\iri2 \ P P 

+ ILd'Vt) + -l-V • (T 2 D"X7x 1 + LVp + AVT) 

1 J an ' 

2 / 2 \ 

+ ( Vittfc + Vfc^ - -5wV- ul V^ttjfe, (86) 

DfMf + p~ x Vip = p _1 V fc ry ^V^u fc + V fc u f - ^fc^V • u^j . (87) 

For the chosen set of fields, n = p/T and p — p [(mi — mz) X\ + 777-2] /T. These 
equations are exact to second order in the spatial gradients for a low density 
Boltzmann gas. Note that in Eqs. (|84|) - ([87|) the second order contributions 
to the cooling rate have been neglected. These second order terms have been 
calculated for a monocomponent fluid [5] and found to be very small relative 
to corresponding terms from the fluxes. Consequently, they have not been 
considered in the hydrodynamic equations (|B3 ]) -([8"7 |) . 

One of the main peculiarities of the granular gases (in contrast to ordinary 
fluids) is the existence of non-trivial solutions to the Navier-Stokes equations 
(|84l) (|87p. even for spatially homogeneous states, 



d t x 1H = 0, d t u H = 0, (88) 

[d t + ((x 1H ,T H ,p H )]T H = 0, [d t + ((x 1H ,T H ,p H )}p H = 0, (89) 

where the subscript H denotes the homogeneous state. Since the dependence 
of the cooling rate C{ x ih,Th,Ph) on xih,Th,Ph is known [Hl|44], these 
first order nonlinear equations can be solved for the time dependence of the 
homogeneous state. The result is the familiar Haff cooling law for T(t) at 
constant density p] HI] : 

As said before, each partial temperature Tj(t) has the same time dependence 
but with a different value [15j . 



where 7 = Tin (£)/T 2 /f (t) is the time-independent temperature ratio. 
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Nevertheless, the homogeneous cooling state (HCS) is unstable to suffi- 
ciently long wavelength perturbations. For systems large enough to support 
such spontaneous fluctuations, the HCS becomes inhomogeneous at long times. 
This feature was first observed in MD simulations of free monocomponent 
gases [55. ■ In MD simulations the inhomogeneities may grow by the forma- 
tion of clusters, ultimately aggregating to a single large cluster [56 ; if cluster 
growth is suppressed, a vortex field may grow to the system size where peri- 
odic boundary conditions can induce a transition to a state with macroscopic 
shear. The mechanism responsible for the growth of inhomogeneities can be 
understood at the level of the Navier-Stokes hydrodynamics, where linear 
stability analysis shows two shear modes and a heat mode to be unstable 



The objective here is to extend this analysis to the case of a binary mixture. 
To do that, we perform a linear stability analysis of the nonlinear hydrody- 
namic equations l[84 p -(|87 p with respect to this HCS for small initial spatial 
perturbations. For ordinary fluids such perturbations decay in time according 
to the hydrodynamic modes of diffusion (shear, thermal, mass) and damped 
sound propagation 3, 58, 59]. For inelastic collisions, the analysis is for fixed 
coefficients of restitution in the long wavelength limit. As will be seen be- 
low, the corresponding modes for a granular mixture are then quite different 
from those for ordinary mixtures. In fact, an alternative study with fixed long 
wavelength and coefficients of restitution approaching unity yields the usual 
ordinary fluid modes. Consequently, the nature of the hydrodynamic modes 
is non uniform with respect to the inelasticity and the wavelength of the 
perturbation. 

Let us assume that the deviations dy a (r,t) = y a (r,t) — yua(t) are small. 
Here, Sy a (r,t) denotes the deviation of {xi,u,T, p} from their values in the 
HCS. If the initial spatial perturbation is sufficiently small, then for some 
initial time interval these deviations will remain small and the hydrodynamic 
equations ([8"4] ) -([8"7 [) can be linearized with respect to Sy a (r,t). This leads to 
a set of partial differential equations with coefficients that are independent 
of space but which depend on time. As in the monocomponent case [57] . 
this time dependence can be eliminated through a change in the time and 
space variables, and a scaling of the hydrodynamic fields. We introduce the 
following dimcnsionlcss space and time variables: 



where vonit) is an effective collision frequency for hard spheres and vqh = 
y2T^r(mi + m-2)lm\m,2- Since {x\h , u# , Tjj ,ph} are evaluated in the HCS, 
then Eqs. (|88|) and (|89|) hold. A set of Fourier transformed dimensionless 
variables are then introduced as 



00 EH [57]. 




s 




(92) 



Sx lk (r) 



Wk(r) 



<5u k (r) 

vqh(t) ' 



ST k (r) 
Th{t) 



Ph(t) ' 
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where Sy^ = {<5xik, Su^, ST^, Sp^} is defined as 



5Vak(r) = J ds 



3 (5?/ Q (s,r). 



(94) 



Note that here the wave vector k is dimcnsionless. 

In terms of the above variables, the transverse velocity components wid = 
wit — (wk • k)k (orthogonal to the wave vector k) decouple from the other 
four modes and hence can be obtained more easily. They obey the equation 



0_ 



i]*k 2 w k 



where £* = Ch/^oh and 



V = - .,2 V, 



0, 



PHV 0H 

where pn = ntinm + m2U2H- The solution for Wkj_(r) reads 



where 



w k ±(r) = w k i(0)exp[s_L(fc)r], 
1 



*;„2 



(95) 
(96) 

(97) 
(98) 



This identifies d — 1 shear (transversal) modes. We see from Eq. (|98|) that 
there exists a critical wave number k ^ given by 



kl = 



2rf 



1/2 



(99) 



This critical value separates two regimes: shear modes with k > k\ always 
decay while those with k <¥]_ grow exponentially. 

The remaining modes are called longitudinal modes. They correspond to 
the set {pk, #k, n^, Wk||} where the longitudinal velocity component (parallel 
to k) is Wk|| = Wk • k. These modes are the solutions of the linear equation 

m 

dSz ak (r) 



dr 



(m<°> + ikM§ + fc 2 Afg) 8z m {r), (100) 

where 6z a k(T) denotes now the four variables {pk, &k, n^, Wk||} ■ The matrices 
in Eq. (|100|) are given by 



M<°> 
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V o o o \cJ 



(101) 
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In these equations, = mx/mz, Xi = riiH/nn, and we have introduced the 
reduced Navier-Stokes transport coefficient^]] 



D* = 



■D, 



n H v 0H 



d ; = 



Ph^oh 



mim 2 n 2 H v 2 H 



D p , D'* = 



Ph^oh 
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n H v. 
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The longitudinal modes have the form exp[s„(fc)r] with n = 1,2,3,4, 
where s n (k) are the eigenvalues of the matrix M(k) = M (0) + ikM^ + fc 2 M (2) , 
namely, they are the solutions of the quartic equation 



det I 



sll =0. 



(109) 



The solution to (|109[) for arbitrary values of k is quite intricate. It is instructive 
to consider first the solutions to these equations in the extreme long wave- 
length limit, k = 0. In this case, they are found to be the eigenvalues of the 
matrix of M< ': 



(110) 



Hence, at asymptotically long wavelengths (k = 0) the spectrum of the lin- 
earized hydrodynamic equations (both transverse and longitudinal) is com- 
prised of a decaying mode at — C*/2, a two-fold degenerate mode at 0, and a 
(i-fold degenerate unstable mode at £*/2. Consequently, some of the solutions 
are unstable. The two zero eigenvalues represent marginal stability solutions, 



1 Note that the definition for the reduced diffusion coefficient D* given here differs 
from the one introduced in Sec. 15.11 
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Fig. 7. Dispersion relations for a = 0.9, xi = 0.2, u = 1 and fi = 4. 

while the negative eigenvalue gives stable solutions. For general initial per- 
turbations all modes are excited. These modes correspond to evolution of 
the fluid due to uniform perturbations of the HCS, i.e., a global change in 
the HCS parameters. The unstable modes are seen to arise from the initial 
perturbations w k j_(0) or iu k ii(0). The marginal modes correspond to changes 
in the composition at fixed pressure, density, and velocity, and to changes in 

— #k at constant composition and velocity. The decaying mode corresponds 
to changes in the temperature or pressure for 77k = #k- The unstable modes 
may appear trivial since they are due entirely to the normalization of the 
fluid velocity by the time dependent thermal velocity. However, this normal- 
ization is required by the scaling of the entire set of equations to obtain time 
independent coefficients. 

The real parts of the modes s±(k) and s n (k) is illustrated in Fig. [7]in the 
case of hard spheres (d = 3) for a = a%j = 0.9, (T\/u2 = 1, sci = 0.2, and 
m\jmi = 4. The k = values correspond to five hydrodynamic modes with 
two different degeneracies. The shear mode degeneracy remains at finite k but 
the other is removed at any finite k. At sufficiently large k a pair of real modes 
become equal and become a complex conjugate pair at all larger wave vectors, 
like a sound mode. The smallest of the unstable modes is that associated with 
the longitudinal velocity, which couples to the scalar hydrodynamic fields. It 
becomes negative at a wave vector smaller than that of Eq. (|99|) and gives the 
threshold for development of spatial instabilities. 

The results obtained here for mixtures show no new surprises relative to 
the case for a monocomponent gas [51 [521 1SZ] , with only the addition of the 
stable mass diffusion mode. Of course, the quantitative features can be quite 
different since there are additional degrees of freedom with the parameter set 
{xih, mi/m.2, Ci/o^j c*ij '}■ Also, the manner in which these linear instabilities 
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are enhanced by the nonlinearities may be different from that for the one 
component case |60j . 

9 Segregation in Granular Binary Mixtures: Thermal 
Diffusion 

The analysis of the linearized hydrodynamic equations for a granular binary 
mixture has shown that the resulting equations exhibit a long wavelength 
instability for d of the modes. These instabilities lead to the spontaneous for- 
mation of velocity vortices and density clusters when the system evolves freely. 
A phenomenon related with the density clustering is the separation or species 
segregation. Segregation and mixing of dissimilar grains is perhaps one of the 
most interesting problems in agitated granular mixtures. In some processes 
it is a desired and useful effect to separate particles of different types, while 
in other situations it is undesired and can be difficult to control. A variety 
of mechanisms have been proposed to describe the separation of particles of 
two sizes in a mixture of vertically shaken particles. Different mechanisms 
include void filling, static compressive force, convection, condensation, ther- 
mal diffusion, interstitial gas forcing, friction, and buoyancy |61| . However, 
in spite of the extensive literature published in the past few years on this 
subject, the problem is not completely understood yet. Among the different 
competing mechanisms, thermal diffusion becomes one of the most relevant 
at large shaking amplitude where the sample of macroscopic grains resembles 
a granular gas. In this regime, binary collisions prevail and kinetic theory can 
be quite useful to analyze the physical mechanisms involved in segregation 
processes. 

Thermal diffusion is caused by the relative motion of the components of a 
mixture because of the presence of a temperature gradient. Due to this motion, 
concentration gradients subsequently appear in the mixture producing diffu- 
sion that tends to oppose those gradients. A steady state is finally achieved 
in which the separation effect arising from thermal diffusion is compensated 
by the diffusion effect. In these conditions, the so-called thermal diffusion fac- 
tor Aij characterizes the amount of segregation parallel to the temperature 
gradient. In this Section, the thermal diffusion factor is determined from the 
Chapman-Enskog solution described before. 

To make some contact with experiments, let us assume that the binary 
granular mixture is in the presence of the gravitational field g = —ge z , where g 
is a positive constant and e z is the unit vector in the positive direction of the z 
axis. In experiments |41j , the energy is usually supplied by vibrating horizontal 
walls so that the system reaches a steady state. Here, instead of considering 
oscillating boundary conditions, particles are assumed to be heated by the 
action of the stochastic driving force ([50]) . which mimics a thermal bath. 
As said above, although the relation between this driven idealized method 
with the use of locally driven wall forces is not completely understood, it 
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must be remarked that in the case of boundary conditions corresponding to 
a sawtooth vibration of one wall the condition to determine the temperature 
ratio coincides with the one derived from the stochastic force [31]. The good 
agreement between theory and simulation found in Fig. [1] for the temperature 
ratio confirms this expectation. 

The thermal diffusion factor Aij (i ^ j) is defined at the steady state in 
which the mass fluxes ji vanish. Under these conditions, the factor A%j is given 
through the relation [62] 

-y^VlnT=— Vx u + Aji = 0. (Ill) 

X i X j 

The physical meaning of Aij can be described by considering a granular binary 
mixture held between plates at different temperatures T (top plate) and T' 
(bottom plate) under gravity. For the sake of concreteness, we will assume 
that gravity and thermal gradient point in parallel directions, i.e., the bottom 
is hotter than the top (T 1 > T). In addition, without loss of generality, we 
also assume that o\ > a 2 - In the steady state, Eq. (|llip describes how the 
thermal field is related to the composition of the mixture. Assuming that A\ 2 is 
constant over the relevant ranges of temperature and composition, integration 
of Eq. (fiTT]) yields 

ln^4 = ^ 12 ln^, (112) 
x 2 x\ T 

where Xi refers to the mole fraction of species i at the top plate and x\ refers 
to the mole fraction of species i at the bottom plate. Consequently, according 
to Eq. (|112p . if A12 > 0, then x' x < Xi, while if A X2 < 0, then x[ > xi. In 
summary, when A12 > 0, the larger particles accumulate at the top of the 
sample (cold plate), while if A X2 < 0, the larger particles accumulate at the 
bottom of the sample (hot plate). The former situation is referred to as the 
Brazil-nut effect (BNE) while the latter is called the reverse Brazil-nut effect 
(RBNE). 

The RBNE was first observed by Hong et al. [63] in MD simulations of 
vertically vibrated systems. They proposed a very simple segregation criterion 
that was later confirmed by Jenkins and Yoon [64] by using kinetic theory. 
More recently, Breu et al. [65] have experimentally investigated conditions 
under which the large particles sink to the bottom and claim that their ex- 
periments confirm the theory of Hong et al. [63] provided a number of condi- 
tions are chosen carefully. In addition to the vertically vibrated systems, some 
works have also focused in the last few years on horizontally driven systems 
showing some similarities to the BNE and its reverse form [66] . However, it 
is important to note that the criterion given in Ref. , 63] is based on some 
drastic assumptions: elastic particles, homogeneous temperature, and energy 
equipartition. These conditions preclude a comparison of the kinetic theory 
derived here with the above simulations. 

Some theoretical attempts to assess the influence of non-equipartition on 
segregation have been recently published. Thus, Trujillo et al. [67] have derived 
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Fig. 8. Phase diagram for BNE/RBNE for mixtures constituted by spheres (d = 3) 
of the same mass density and equal total volumes of large and small particles. The 
data points represent the MD simulation results [41] for a = 0.78 when convection 
is suppressed. Points below (above) the curve correspond to RBNE (BNE). 

an evolution equation for the relative velocity of the intruders starting from 
the kinetic theory proposed by Jenkins and Yoon [64 , which applies for weak 
dissipation. They use constitutive relations for partial pressures that take into 
account the breakdown of energy equipartition between the two species. How- 
ever, the influence of temperature gradients, which exist in the vibro-fluidized 
regime, is neglected in Ref. [67] because it is assumed that the pressure and 
temperature are constant in the absence of the intruder. A more refined theory 
has recently been provided by Brey et al. [42] in the case of a single intruder 
in a vibrated granular mixture under gravity. The theory displayed in this 
section covers some of the aspects not accounted for in the previous theories 
[42] [64] [67] since it is based on a kinetic theory [19] that goes beyond the 
quasi-elastic limit [64] [67] and applies for arbitrary composition x\ (and so, 
it reduces to the results obtained in Ref. [15] when x\ — > 0). This allows one 
to assess the influence of composition and dissipation on thermal diffusion in 
bi-disperse granular gases without any restriction on the parameter space of 
the system. 

To determine the dependence of the coefficient A\2 on the parameters of 
the mixture, we consider a non-convecting (u = 0) steady state with only 
gradients along the vertical direction (z axis). In this case, the mass balance 
equation (fl~2"]) yields ji = j 2 = 0, while the momentum equation (fl~3"]) gives 



dp 

dz 



-pg- 



(113) 



To first order in the spatial gradients, the constitutive equation for the mass 
flux ji jZ is given by Eq. ([55]) . i.e., 
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where the susceptibility coefficient Xij = in the particular case of the grav- 
itational force. The condition j\ z — yields 
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dz m\minp D 



D' dT 



m,\mip D dz ' 



(115) 



where use has been made of Eq. (|113p . Substitution of Eq. (|115[) into Eq 
(TTTTj) leads to 

np 2 D' - D p g* 



where 



An = =r^, (116) 

P1P2 D 

n {01 /az) 

is the reduced gravity acceleration. Since the mutual diffusion coefficient D 
is positive [HI 03], * ne sl S n °f ^12 is determined by the sign of the quantity 
D' — D p g* . This result is general since it goes beyond the regime of density 
considered. 
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Fig. 9. Phase diagram for BNE/RBNE in three dimensions for Qy = 0.7 and three 
values of composition: (a) x\ — 0, (b) xi = 0.3, and (c) x\ — 0.7. Points below 
(above) each curve correspond to RBNE (BNE). 



To gain some insight into the explicit dependence of D' and D p on the 
parameter space of the system, one has to resort to a kinetic theory descrip- 
tion. For a low-density gas, the expressions of the coefficients D 1 and D p in 
the first-Sonine approximation are given by 
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Fig. 10. Phase diagram for BNE/RBNE in three dimensions for xi = i and three 
values of the (common) coefficient of restitution: (a) a = 0.9, (b) a — 0.8, and (c) 
a — 0.5. Points below (above) each curve correspond to RBNE (BNE). 



D' = 0, D p 



PiP x 2 I 7 

p 2 V X 2 + Xi~f \/J, 



where v is the (positive) collision frequency [44) 

2^ (d - 1)/2 d -i n 
"- — m 12 w (l + a 12 ) 



1/2 



(118) 



(119) 



Given that the driving stochastic term does not play a neutral role in the 
transport, it must be remarked that the expressions for the transport coeffi- 
cients obtained in the driven case slightly differ from the ones derived in the 
free cooling case [THl HI] . 

Consequently, according to Eqs. and (|118[) . the sign of A\2 is the 

same as that of the pressure diffusion coefficient D p . The condition A12 = 
(or equivalently, D p = 0) provides the criterion for the transition from BNE 
to RBNE. Equation (|118[) shows that the sign of D p is determined by the 
value of the control parameter 



2 



miT 2 



(120) 



This parameter gives the mean square velocity of the large particles relative 
to that of the small particles. Thus, when 9 > 1 (9 < 1), the thermal diffusion 
factor is positive (negative), which leads to BNE (RBNE). The criterion for 
the transition condition from BNE to RBNE is 9 = 1, i.e., 



mi 
m 2 



T 2 



(121) 
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In the case of equal granular temperatures (energy equipartition) , 9 — > 
and so segregation is predicted for particles that differ in mass, no matter 
what their diameters may be [53]. It must be remarked that, due to the lack 
of energy equipartition, the condition 6 = 1 is rather complicated since it 
involves all the parameter space of the system. In particular, even when the 
species differ only by their respective coefficients of restitution they also seg- 
regate when subject to a temperature gradient. This is a novel pure effect of 
inelasticity on segregation [25 , 68J. On the other hand, the criterion (|121|) for 
the transition BNE<=>RBNE is the same as the one found previously in Ref. 
[67] when otij is close to 1 and in Ref. [42] in the intruder limit case (x\ — > 0). 
However, as said before, the results obtained here are more general since they 
cover all the range of the parameter space of the system. 

To illustrate size segregation driven by thermal diffusion, we consider mix- 
tures constituted by spheres (d = 3) of the same material and equal total 
volumes of large and small particles. In this case, m\jm2 — (C1/CT2) 3 and 
xz/xi = (ci/o^) 3 . Figure [8] shows the phase diagram BNE/RBNE for this 
kind of systems. The data points represent the simulation results obtained by 
Schroter et al. [UJ for a = 0.78 in agitated mixtures constituted by parti- 
cles of the same density. To the best of my knowledge, this is one of the few 
experiments in which thermal diffusion has been isolated from the remain- 
ing segregation mechanisms [61]. Our results show that, for a given value of 
the coefficient of restitution, the RBNE is dominant at small diameter ra- 
tios. However, since non-equipartition grows with increasing diameter ratio, 
the system shows a crossover to BNE at sufficiently large diameter ratios. 
This behavior agrees qualitatively well with the results reported in Ref. [4T] 
at large shaking amplitudes, where thermal diffusion becomes the relevant 
segregation mechanism. At a quantitative level, we observe that the results 
are also consistent with the simulation results reported in [UJ when periodic 
boundary conditions are used to suppress convection since they do not ob- 
serve a change back to BNE for diameter ratios up to 3 (see red squares in 
Fig. 11 of [UJ). Although the parameter range explored in MD simulations is 
smaller than the one analyzed here, one is tempted to extrapolate the simula- 
tion data presented in Ref. [41 to roughly predict the transition value of the 
diameter ratio at a = 0.78 (which is the value of the coefficient of restitution 
considered in the simulations). Thus, if one extrapolates from the simulation 
data at the diameter ratios of <J\joi = 2 and <J\/u2 — 3, one sees that the 
transition from RBNE to BNE might be around <ti/(T2 = 10, which would 
quantitatively agree with the results reported in Fig. [8] Figure [8] also shows 
that the BNE is completely destroyed in the quasielastic limit (a ~ 1). 

Let us now investigate the influence of composition on segregation. Figured] 
shows a typical phase diagram in the three-dimensional case for a»j = a = 0.7 
and three different values of the mole fraction x\. The lines separate the 
regimes between BNE and RBNE. We observe that the composition of the 
mixture has significant effects in reducing the BNE as the concentration of 
larger particles increases. In addition, for a given value of composition, the 
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transition from BNE to RBNE may occur following two paths: (i) along the 
constant mass ratio m 2 /mi with increasing size ratio o^/o^, and (ii) along 
the constant size ratio with increasing mass ratio m 2 / 'mi. The influence of 
dissipation on the phase diagrams BNE/RBNE is illustrated in Fig. [TU1 for 
d = 3 in the case of an equimolar mixture (x% = i) and three values of 
the (common) coefficient of restitution a. We observe that the role played by 
inelasticity is quite important since the regime of RBNE increases significantly 
with dissipation. Similar results are found for other values of composition. 

In summary, thermal diffusion (which is the relevant segregation mech- 
anism in agitated granular mixtures at large shaking amplitudes) can been 
analyzed by the Boltzmann kinetic theory. This theory is able to explain 
some of the experimental and/or MD segregation results [41] observed within 
the range of parameter space explored. A more quantitative comparison in 
the dilute regime with MD simulations is needed to show the relevance of 
the Boltzmann equation to analyze segregation driven by a thermal gradi- 
ent. As said before, comparison with MD simulations in the tracer limit case 
(xx — > 0) [32] for a dilute gas has shown the reliability of the inelastic Boltz- 
mann equation to describe segregation. In this context, one expects that the 
same agreement observed before in the intruder case |42j is maintained when 
X\ is different from zero. 



10 Steady States: Uniform Shear Flow 

In the preceding sections, the Navier-Stokes equations (constitutive equa- 
tions that are linear in the hydrodynamic gradients) have been shown to be 
quite useful to describe appropriately several problems in granular mixtures. 
However, under some circumstances large gradients occur and more complex 
constitutive equations are required. The need for more complex constitutive 
equations does not signal a breakdown of hydrodynamics 69J , only a failure 
of the Navier-Stokes approximation 70J . Although in this case the Chapman- 
Enskog method can be carried out to second order in gradients (Burnett or- 
der), it is likely that failure of the Navier-Stokes description signals the need 
for other methods to construct the normal solution that are not based on a 
small gradient expansion. 

One of the most interesting problems in granular fluids is the simple or 
uniform shear flow (USF) [T31[25]. As described in Sec. 03 this state is charac- 
terized by uniform density and temperature and a simple shear with the local 
velocity field given by u\ tX = U2 tX = ay, Uy = u z — 0, where a is the constant 
shear rate. The USF is a well-known nonequilibrium problem widely studied, 
for both granular [26l EH ED EH [73] and conventional [28] [74] gases. However, 
the nature of this state is quite different in both systems since a steady state 
is achieved for granular fluids when viscous (shear) heating is compensated 
for by energy dissipation in collisions: 
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(122) 



This steady state is what we want to analyze in this section. The balance 
equation (|122[) shows the intrinsic connection between the shear field and dis- 
sipation in the system. This contrasts with the description of USF for elastic 
fluids where a steady state is not possible unless an external thermostat is in- 
troduced [2E] ■ Note that the hydrodynamic steady shear flow state associated 
with the condition (|122|) is inherently beyond the scope of the Navier-Stokes 
or Newtonian hydrodynamic equations 29J. The reason for this is the exis- 
tence of an internal mechanism, collisional cooling, that sets the strength of 
the velocity gradient in the steady state. For normal fluids, this scale is set by 
external sources (boundary conditions, driving forces) that can be controlled 
to admit the conditions required for Navier-Stokes hydrodynamics. In con- 
trast, collisional cooling is fixed by the mechanical properties of the particles 
making up the fluid. This observation is significant because it prevents the 
possibility of measuring the Newtonian shear viscosity for granular fluids in 
the steady USF [721 E3] • More generally, it provides a caution regarding the 
simulation of other steady states to study Navier-Stokes hydrodynamics when 
the gradients are strongly correlated to the collisional cooling [29 . 

From a microscopic point of view, the simple shear flow problem becomes 
spatially uniform in the local Lagrangian frame moving with the flow velocity 
u. In this frame [2H1 HS1 E3) the velocity distribution functions adopt the 
form: /j(r, v) — > /i(V), where Vk = Vf, — a^t is the peculiar velocity. Here, 
o>M = adkxfiiy Under these conditions, the set of Boltzmann kinetic equations 
(with T{ = 0) for an isolated system reads 



The most relevant transport properties in a shear flow problem are obtained 
from the pressure tensor P = Pi + P2, where P; is the partial pressure tensor 
of the species i given by 



The trace of Pj defines the partial temperatures Ti as Tj = TrP i/drii. 
As said before, these temperatures measure the mean kinetic energy of each 
species. The elements of the pressure tensor Pj can be obtained by multiplying 
the Boltzmann equation (| 1 23[) by m^VV and integrating over V. The result 



aV v7 —m) = £ J tf [V|/i, fj] , (i = 1, 2). 
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where we have introduced, the collisional moments a-s 
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(126) 



From Eq. (|125p . in particular, one gets the balance equation for the partial 
temperature T, 

aPi, xy = --jPidi (127) 

where pi = n{Ti is the partial pressure of species i and Q is defined by Eq. 
(fT0|) . According to Eq. (|127p . the (steady) partial temperature in the simple 
shear flow problem can be obtained by equating the viscous heating term 
a\Pi,xy \ to the collisional cooling term (d/2)piQ. 

The determination of requires the knowledge of the velocity distribu- 
tion functions /j. This is quite a formidable task, even in the monocomponcnt 
case [27] ■ However, as in the elastic case, one expects to get a good estimate 
of Aij by using Grad's approximation [15] : 



/(v) -> f iiM (y) ( i + ^±c i>M v k Vi 



(128) 



where fi_M is a Maxwellian distribution at the temperature of the species i, 
i.e., 

rrii \ ' I miV 



fi,M(V) = rii 



\2irTi 



exp 



221 



(129) 



As happens in the case of homogeneous states, in general the three temper- 
atures T, Ti, and are different in the inelastic case. For this reason we 
choose the parameters in the Maxwellians so that it is normalized to rii and 
provides the exact second moment of fi. The Maxwellians /^m for the two 
species can be quite different due to the temperature differences. This aspect 
is essential in a two-temperature theory and has not been taken into account 
in most of the previous studies (TU [711 [72] . The coefficient Q can be identified 
by requiring the moments with respect to VV of the trial function (|128p to 
be the same as those for the exact distribution fi. This leads to 



Pi 
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With this approximation, the Boltzmann collisional moments A^- can be 
explicitly evaluated. The result is [TTJ [78] 
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where 
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Tj - Ti 



HJ " (mj/m^Ti + Tj dT3 (M + 3 ~ 3 " ij ; ' 
The partial cooling rates can be easily obtained from Eqs. (fT0|) and (] 1 3 1 [) 



(132) 
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Fig. 11. Plot of the temperature ratio T1/T2 as a function of the size ratio 0-1/02 = 
(mi/m^) 1 / 2 for a two-dimensional system in the case x\ = 1/2 and two different 
values of the (common) coefficient of restitution: a — 0.9 and a = 0.8. The solid 
lines are the theoretical predictions based on Grad's solution, while the symbols 
refer to the DSMC results. The dashed lines correspond to the results obtained from 
the stochastic thermostat condition (1311). 



Substitution of Eq. (|131[) into the set of equations (|125[) allows one to get 
the partial pressure tensor P, in terms of the temperature ratio 7 = T1/T2 
and the parameters of the mixture. The temperature ratio can be obtained 
from Eq. (fT27|> as 

X2(,2Pl,xy 



7 



(133) 
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When the expressions of Pi and Q are used in Eq. Q133p , one gets a closed 
equation for the temperature ratio 7, that can be solved numerically. In Fig. 
1111 we plot 7 versus the diameter ratio o\ / o~2 for a two-dimensional (d — 
2) granular gas with x\ = 1/2 and two different values of a. The symbols 
refer to the simulation data obtained from the DSMC method [75]. Here, 
we have assumed that the disks are made of the same material, and hence 
oiij = a and 7711/7712 = (ci/o^) 2 . The dependence of 7 on eri/<72 obtained in 
the homogeneous steady state driven by the stochastic thermostat (|3Tj|) is also 
included for comparison. It is clearly seen that the kinetic theory results based 
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on Grad's solution agree very well with simulation data, even for quite large 
values of the size ratio. In addition, the thermostat results overestimate the 
simulation ones (especially for large mass ratio), showing that the properties 
of the system are not insensitive to the way at which the granular gas is driven. 
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Fig. 12. Plot of the reduced pressure p* versus the mass ratio /i = mxjm-x for a two- 
dimensional system with o\ = 02, x\ = 1/2 and a — 0.9. The solid line corresponds 
to the theoretical predictions derived from Grad's solution, the dotted line refers 
to the latter theory but using the expression of T1/T2 obtained from the stochastic 
thermostat condition (|31[) . and the dashed line is the result obtained from Grad's 
solution by assuming the equality of the partial temperatures (7 = 1). The symbols 
are the DSMC results. 




Let us now consider the transport coefficients. To analyze the Theologi- 
cal properties in the steady state, it is convenient to introduce dimensionless 
quantities. As usual [75], for a low-density gas we introduce the reduced pres- 
sure p* and the reduced shear viscosity rf as 

P* = - E ^a> (134) 

T) = — , (135) 

pivfia 



where 77 = —P xy /a is the non-Newtonian shear viscosity, P xy — P\, xy + Pi^xy 
and v = [7r( d - 1 )/ 2 / J T(d/2)]ncr^ 1 u . In Figs. [TJ and HI we plot p* and rf , 
respectively, as functions of the mass ratio fj, = m\jmi for an equal-size 
(or = a 2) binary mixture of disks (d = 2) with x\ — 1/2 and a = 0.9. 



Kinetic Theory for Binary Granular Mixtures at Low-Density 



41 



2 



.00 



T 1 1 1 — I I | 



T — IT 




1 



10 



100 



Fig. 13. Plot of the reduced shear viscosity n* versus the mass ratio fj, = mi/ffl2 
for a two-dimensional system with cti = (72, Si = 1/2 and a = 0.9. The solid line 
corresponds to the theoretical predictions derived from Grad's solution, the dotted 
line refers to the latter theory but using the expression of Ti/T? obtained from the 
stochastic thermostat condition (|31|) . and the dashed line is the result obtained from 
Grad's solution by assuming the equality of the partial temperatures (7 = 1). The 
symbols are the DSMC results. 

We have also included the predictions for p* and 77* given by the kinetic 
theory but taking the expression of 7 derived when the system is driven by 
the stochastic thermostat [ID]. We observe again in both figures an excellent 
agreement between the Boltzmann theory based on Grad's solution and the 
DSMC results, even for very disparate values of the mass ratio. With respect 
to the influence of energy nonequipartition, Fig. [12] shows that p* presents a 
non-monotonic behaviour with the mass ratio whereas the theoretical predic- 
tions with the equipartition assumption monotonically increase with /j. In the 
case of the shear viscosity, as seen in Fig. [T31 both theories (with and with- 
out energy nonequipartition) predict a non-monotonic dependence of rj* on 
/z. However, at a quantitative level, the influence of energy nonequipartition 
is quite significant over the whole range of mass ratios considered. The non- 
monotonic dependence of p* and rj* on \x obtained here from the Boltzmann 
kinetic theory also agrees qualitatively well with MD simulations carried out 
for bidisperse dense systems |72j . Thus, for instance, the minimum values of p* 
and 77* are located close to (1 — 10 in both dilute and dense cases. Moreover, 
the predictions for the transport properties given from the present theory 
by taking the stochastic thermostat expression of 7 are quite close to those 
obtained from the actual value of 7, especially for large mass ratios. 
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11 Summary and Concluding Remarks 

The primary objective of this review has been to derive the Navier-Stokes hy- 
drodynamic equations of a binary mixture of granular gases from the (inelas- 
tic) Boltzmann kinetic theory. The Chapman-Enskog method [4j [43] is used 
to solve the Boltzmann equation up to the first order in the spatial gradients 
and the associated transport coefficients are given in terms of the solutions of 
a set of linear integral equations. These equations have been approximately 
solved by taking the leading terms in a Sonine polynomial expansion. Compar- 
ison with controlled numerical simulations in some idealized conditions shows 
quite a good agreement between theory and simulation even for strong dis- 
sipation. This supports the idea that the hydrodynamic description (derived 
from kinetic theory) appears to be a powerful tool for analysis and predictions 
of rapid flow gas dynamics of polydisperse systems [T2] . 

The reference state in the Chapman-Enskog expansion has been taken 
to be an exact solution of the uniform Boltzmann equation. An interesting 
and important result of this solution [15] is that the partial temperatures 
(which measure the mean kinetic energy of each species) are different. This 
does not mean that there are additional degrees of freedom since the partial 
temperatures can be expressed in terms of the global temperature. This is 
confirmed by noting that Haff 's cooling law [1] (in the free cooling case) is the 
hydrodynamic mode at long wavelengths and MD simulations confirm that 
the global temperature dominates after a transient period of a few collision 
times 34 . In this case, only the global temperature should appear among 
the hydrodynamic fields. Nevertheless, the species temperatures play a new 
and interesting secondary role [20] . For an ordinary (molecular) gas, there is 
a rapid velocity relaxation in each fluid cell to a local equilibrium state on 
the time scale of a few collisions (e.g., as illustrated by the approach to Haft's 
law). Subsequently, the equilibration among cells occurs via the hydrodynamic 
equations. In each cell the species velocity distributions are characterized by 
the species temperatures. These are approximately the same due to equiparti- 
tion, and the hydrodynamic relaxation occurs for the single common temper- 
ature [43] . A similar rapid velocity relaxation occurs for granular gases in each 
small cell, but to a universal state different from local equilibrium and one for 
which equipartition no longer occurs. Hence, the species temperatures Ti are 
different from each other and from the overall temperature T of the cell. Nev- 
ertheless, the time dependence of all temperatures (in the free cooling case) 
is the same in this and subsequent states, i.e., they are proportional to the 
global temperature. This implies that the species temperatures do not provide 
any new dynamical degree of freedom at the hydrodynamic stage. However, 
they still characterize the shape of the partial velocity distributions and affect 
the quantitative averages calculated with these distributions. The transport 
coefficients for granular mixtures therefore have new quantitative effects aris- 
ing from the time independent temperature ratios for each species |19j. This 
view contrasts with some recent works [80 , where additional equations for 
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each species temperature have been included among the hydrodynamic set. 
However, as mentioned before, this is an unnecessary complication, describing 
additional kinetics beyond hydrodynamics that is relevant only on the time 
scale of a few collisions. 

Another important issue discussed here has been the applicability of the 
Navier-Stokes transport coefficients since their expressions are not restricted 
to weak inelasticity [12] . However, the Navier-Stokes hydrodynamic equations 
themselves may or may not be limited with respect to inelasticity, depending 
on the particular states studied. The Chapman-Enskog method assumes that 
the relative changes of the hydrodynamic fields over distances of the order 
of the mean free path are small. In the case of ordinary fluids this can be 
controlled by the initial or boundary conditions. For granular gases the situ- 
ation is more complicated since in some cases (e.g., steady states such as the 
simple shear flow problem [2 9) ) the boundary conditions imply a relationship 
between dissipation and gradients so that both cannot be chosen indepen- 
dently. In these cases, the Navier-Stokes approximation only holds for nearly 
elastic particles. However, the transport coefficients characterizing the Navier- 
Stokes hydrodynamic equations are nonlinear functions of the coefficients of 
restitution, regardless the applicability of those equations. 

In spite of the above cautions, the Navier-Stokes approximation is ap- 
propriate and accurate for a wide class of flows. One group refers to spatial 
perturbations of the homogeneous cooling state (HCS) for an isolated sys- 
tem. Both MD and DSMC simulations [7] have confirmed the dependence 
of the Navier-Stokes transport coefficients on the coefficient of restitution, 
and application of the Navier-Stokes hydrodynamics with these coefficients 
to describe cluster formation has also been confirmed quantitatively [60 . The 
same kinetic theory results apply to driven systems as well. This is so since 
the reference state is a local HCS whose parameters vary throughout the sys- 
tem to match the physical values in each cell. Examples include application of 
Navier-Stokes hydrodynamics from kinetic theory to symmetry breaking and 
density/temperature profiles in vertical vibrated gases, for comparison with 
simulation [81] . Similar comparison with Navier-Stokes hydrodynamics of the 
latter and of supersonic flow past a wedge in real experiments has been given 
|82[ 183] . showing both qualitative and quantitative agreement. In summary, 
the Navier-Stokes equations with the constitutive equations presented here 
remain an important and useful description for a wide class of granular flow, 
although more limited than for normal gases. 

The explicit knowledge of the transport coefficients and the cooling rate 
allows one to make some applications of the Navier-Stokes hydrodynamic 
equations. One of them has been to obtain the linear hydrodynamic equa- 
tions for small perturbations of the homogenous cooling state. The resulting 
equations exhibit a long wavelength instability for three of the modes. This 
is quite similar to the case of a monocomponent granular gas [5j [SH [57] , and 
in fact the same modes are unstable here. The additional diffusion mode for 
two species behaves as for a normal fluid. 
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On the other hand, the constitutive equations for the mass and heat fluxes 
of a granular binary mixture differ from those obtained for ordinary fluids [I] . 
This is because the usual restrictions of irreversible thermodynamics no longer 
apply. These restrictions include Onsager's reciprocal relations among various 
transport coefficients and the extent to which these are violated has also been 
shown here. Another application of the Navier-Stokes equations has been to 
assess the violation of the Einstein relation between the diffusion and mobility 
coefficients. In the undriven case, the analysis shows that this violation is 
due to three independent reasons [53]: the absence of the Gibbs state, the 
cooling of the reference state, and the occurrence of different temperatures for 
the particle and surrounding fluid. However, when the mixture is subjected 
to stochastic driving, a modified Einstein relation suggested by recent MD 
simulations [52] has also been analyzed. In this case, the results show that 
the deviations of the (modified) Einstein ratio from unity are in general very 
small (less than 1%), in agreement with MD simulations 52J. 

Thermal diffusion becomes the relevant segregation mechanism in agitated 
granular mixtures at large shaking amplitudes. In these conditions, the use of 
the Boltzmann kinetic theory for low-density gases appears justified to under- 
stand the influence of thermal gradient on segregation phenomena. The ther- 
mal diffusion factor in a heated granular mixture has been explicitly evaluated 
from the Chapman-Enskog solution to the Boltzmann equation. The results 
show that the criterion for the transition Brazil-nut effect -<=>reverse Brazil- 
nut effect is provided by the control parameter 9 = 77i 2 Ti/miT 2 [25 ] 142 ] l67|. 
Given that the energy equipartition is broken, the condition 9 = 1 is quite 
complex since it involves all the parameters of the system: composition, 
masses, sizes, and coefficients of restitution. The Boltzmann kinetic theory 
results agree qualitatively well with recent MD simulations [41] within the 
range of parameter space analyzed. 

The hydrodynamic description also seems to be justified in the case of 
steady states that are inherently beyond the scope of the Navier-Stokes hy- 
drodynamic equations. The reason for this non-Newtonian behavior is the 
existence of an internal mechanism, collisional cooling, that sets the scale of 
the spatial gradients in the steady state. For ordinary fluids, this scale can be 
externally controlled by external sources so that the conditions for Navier- 
Stokes hydrodynamics apply. On the other hand, for granular gases, collisional 
cooling is fixed by the mechanical properties of the particles of the system and 
so the gas can depart from the Navier-Stokes description. One well-known ex- 
ample of steady states is the simple or uniform shear flow (USF) . However, in 
spite of the extensive prior work on USF for granular fluids [271 E2 E3] , the in- 
herent non-Newtonian character of this state has not been conveniently taken 
into account. In fact, MD simulations of steady USF have been used for gran- 
ular fluids to measure the Newtonian or Navier-Stokes shear viscosity. The 
results derived here from Grad's solution and DSMC simulations show that 
USF is an ideal testing ground for the study of rheology since any choice of the 
shear rate and the coefficients of restitution will provide non-Newtonian 
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effects. It is one of the fascinating features of granular fluids that phenom- 
ena associated with complex fluids are more easily accessible than for simple 
atomic fluids [T21I53] . 

Hydrodynamics derived from hard-sphere models have found widespread 
use in the description of numerous industrial processes involving solid parti- 
cles. Of particular relevance are high-speed, gas-solid flows as found in pneu- 
matic conveyors (of ores, chemicals, grains, etc.) and fluidized beds (for fluid 
catalytic cracking, power generation, granulation of pharmaceutical powders, 
synthesis of fine chemicals like titania, etc.). Such descriptions are now stan- 
dard features of commercial and research codes. Those codes rely upon accu- 
rate transport properties and a first order objective is to assure this accuracy 
from a careful theoretical treatment. As shown in this review, the price of this 
approach, in contrast to more phcnomcnological approaches, is an increasing 
complexity of the expressions as the systems become more complex. 

The analysis carried out in this presentation has been focused on mix- 
tures in the dilute regime, where the collisional transfer contributions to the 
transport coefficients are neglected and only their kinetic contributions are 
considered. A further step is to develop a theory for moderately dense gran- 
ular mixtures. This will provide a fundamental basis for the application of 
hydrodynamics under realistic conditions. Possible extension of the present 
Boltzmann kinetic theory to higher densities can be done in the context of 
the Enskog kinetic equation [43] . Preliminary results [85] have been restricted 
to the uniform shear flow state to get directly the shear viscosity coefficient of 
a heated granular mixture. The extension of this study [85] to states with gra- 
dients of concentration, pressure, and temperature is somewhat intricate due 
to subtleties associated with the spatial dependence of the pair correlations 
functions considered in the revised Enskog theory. A future work is to ex- 
tend the results derived for moderately dense mixtures of smooth elastic hard 
spheres [86] to inelastic collisions. This would allow us to assess the influence 
of density on the different problems addressed in this review. Of course, the 
precise expressions for transport coefficients in this case will be even more 
complex than for a dilute gas due to the expanded parameter space. However, 
this complexity is not a problem for implementation in a code. 

As shown along this overview, granular mixtures exhibit a wide range of 
interesting phenomena for which the Navier-Stokes hydrodynamic equations 
can be considered as an accurate and practical tool. However, due to their 
complexity, many of their features are not fully understood. Kinetic theory 
and hydrodynamics (in the broader sense) can be expected to provide some 
insight into the understanding of such complex materials. 
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